clear;
clc;
addpath(genpath('../util'))
%大尺度衰弱参数
D0=1;
C0=10^(-30/10);
exponent_BS_LIS=2.2;
exponent_LIS_USER=2.8;

%BS位置
x_BS=0;
y_BS=0;
z_BS=25;

%LIS位置
x_IRS=0;
y_IRS=10;
z_IRS=40;

%用户位置范围
x_left=-50;
x_right=50;
y_left=10;
y_right=110;

%噪声功率
sigma=10^(-80/10)*10^(-3);

%BS-LIS距离
d_BS_LIS=sqrt((x_BS-x_IRS)^2+(y_BS-y_IRS)^2);
[pathGain_BS_LIS] = pathLOS(d_BS_LIS,D0,exponent_BS_LIS,C0);

L = 0;
K = 4;
M=8;
F=eye(L);
NUM_timesamples = 1;
pathGain = ones(K,1);
K_factor = 0*ones(K,1);
e=10^(-3);
%disp('H1');
for ueIdx = 1:K
    tau(ueIdx)= 0; %rand;sqrt(0.01)
end
p_correlation = 0.3*ones(K,1);
for k = 1:K  
    for kk=1:L
        for ll=1:L
            T2_base(kk,ll,k) = (p_correlation(k))^(abs(kk-ll));
        end
    end
end


for k=1:K
    x_user=rand*(x_right-x_left)+x_left;
    y_user=rand*(y_right-y_left)+y_left;
    d_IRS_user=sqrt((x_user-x_IRS)^2+(y_user-y_IRS)^2);
    pathGain(k,1)=pathLOS(d_IRS_user,D0,exponent_LIS_USER,C0);
end
H1_base(:,:,1)=randn(L,M)+j*randn(L,M);
H1(:,:,1) = H1_base(:,:,1)*sqrt(L*pathGain_BS_LIS/trace(H1_base(:,:,1)*H1_base(:,:,1)'));
H1(:,:,2) = H1_base(:,:,1)*sqrt(L*pathGain_BS_LIS*10^(7)/trace(H1_base(:,:,1)*H1_base(:,:,1)'));
Hbar2_base(:,:,1) =randn(K,L)+j*randn(K,L);
for ueIdx = 1:K
    T2(1:L,1:L,ueIdx,1) = pathGain(ueIdx)*(L/trace(T2_base(1:L,1:L,ueIdx)))*T2_base(1:L,1:L,ueIdx)/(K_factor(ueIdx)+1);
    Hbar2(ueIdx,1:L,1) = Hbar2_base(ueIdx,1:L,1) * sqrt(pathGain(ueIdx)) * sqrt(1/trace(Hbar2_base(ueIdx,1:L,1)*Hbar2_base(ueIdx,1:L,1)')) * sqrt(K_factor(ueIdx)/(K_factor(ueIdx)+1)) ;
end

for ueIdx = 1:K
    T2(1:L,1:L,ueIdx,2) = pathGain(ueIdx)*10^(7)*(L/trace(T2_base(1:L,1:L,ueIdx)))*T2_base(1:L,1:L,ueIdx)/(K_factor(ueIdx)+1);
    Hbar2(ueIdx,1:L,2) = Hbar2_base(ueIdx,1:L,1) * sqrt(pathGain(ueIdx)*10^(7)) * sqrt(1/trace(Hbar2_base(ueIdx,1:L,1)*Hbar2_base(ueIdx,1:L,1)')) * sqrt(K_factor(ueIdx)/(K_factor(ueIdx)+1)) ;
end


almax = 10^4;
almin = 10^(-8);

alpha_up = almax; %/rho
alpha_low = almin;%/rho

F_temp=eye(L);
P=10;
d=5;
f=5;
c1=2;
c2=2;
NUM=70;
result_PSO1=zeros(1,(P-f)/d+1);
result_PSO2=zeros(1,(P-f)/d+1);
result_PSO3=zeros(1,(P-f)/d+1);

%发送功率
power1=10^(45/10)*10^(-3);


for i=1:NUM_timesamples
    i
    for p = f:d:P
        p
        rho1 = power1/sigma;
       
       [ F_best_PSO1,max_PSO1,alpha1 ]=AO(p,c1,c2,L,NUM,K,T2(:,:,:,1),Hbar2(:,:,1),tau,rho1,M,H1(:,:,1),e,"fix",alpha_low,alpha_up);
       [ F_best_PSO2,max_PSO2,alpha2 ]=AO(p,c1,c2,L,NUM,K,T2(:,:,:,2),Hbar2(:,:,2),tau,rho1/10^(14),M,H1(:,:,2),e,"fix",alpha_low,alpha_up);
       max_PSO1
       max_PSO2
       
       Sumrate1 = Sumrate_MonteCarlo_test(K,L,M,rho1,1/rho1,1000,T2(:,:,:,1),Hbar2(:,:,1),eye(L),H1(:,:,1),tau)
       Sumrate2 = Sumrate_MonteCarlo_test(K,L,M,rho1/10^(14),1/rho1*10^(14),1000,T2(:,:,:,2),Hbar2(:,:,2),eye(L),H1(:,:,2),tau)

       
    end
end
